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Abstract 

We present a new model and a novel loosely coupled partitioned numeri- 
cal scheme modeling fluid-structure interaction (FSI) in blood flow allowing 
non-zero longitudinal displacement. Arterial walls are modeled by a linearly 
viscoelastic, cylindrical Koiter shell model capturing both radial and longi- 
tudinal displacement. Fluid flow is modeled by the Navier-Stokes equations 
for an incompressible, viscous fluid. The two are fully coupled via kine- 
matic and dynamic coupling conditions. Our numerical scheme is based on 
a new modified Lie operator splitting that decouples the fluid and structure 
sub-problems in a way that leads to a loosely coupled scheme which is uncon- 
ditionally stable. This was achieved by a clever use of the kinematic coupling 
condition at the fluid and structure sub-problems, leading to an implicit cou- 
pling between the fluid and structure velocities. The proposed scheme is a 
modification of the recently introduced "kinematically coupled scheme" for 
which the newly proposed modified Lie splitting significantly increases the 
accuracy. The performance and accuracy of the scheme were studied on 
a couple of instructive examples including a comparison with a monolithic 
scheme. It was shown that the accuracy of our scheme was comparable to 
that of the monolithic scheme, while our scheme retains all the main ad- 
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vantages of partitioned schemes, such as modularity, simple implementation, 
and low computational costs. 

Keywords: Fluid-structure interaction, loosely coupled scheme, blood flow. 



1. Introduction 

We study fluid-structure interaction (FSI) between an incompressible vis- 
cous, Newtonian fluid, and a thin viscoelastic structure modeled by the lin- 
early viscoelastic cylindrical Koiter shell model. The cylindrical viscoelastic 
Koiter shell model is derived to describe the mechanical properties of arterial 
walls, while the Navier-Stokes equations for an incompressible, viscous, New- 
tonian fluid were employed to model the flow of blood in medium-to-large 
human arteries. The two are coupled via the kinematic (no-slip) and dynamic 
(balance of contact forces) coupling conditions. Motivated by recent results 
of in vivo measurements of arterial wall motion [H El EJ H] , which indicate 
that both the radial and longitudinal displacement, as well as viscoelasticity 
of arterial walls, are important in disease formation, we derived in this work 
the viscoelastic cylindrical Koiter shell model which captures both radial and 
longitudinal displacement, with the viscoelasticity of Kelvin- Voigt type. The 
novel Koiter shell model is then coupled to the Navier-Stokes equations, and 
the coupled FSI problem is solved numerically. In this manuscript we devise 
a stable, loosely coupled scheme to numerically solve the fully coupled FSI 
problem. The scheme is based on a novel modified Lie's time-splitting, and on 
an implicit use of the kinematic coupling condition, as in [5], which provides 
stability of the scheme without the need for sub-iterations between the fluid 
and structure sub-solvers. Stability of the scheme was proved in [6] on the 
same, simplified benchmark problem as in [TJ. We provide numerical results 
which show that the scheme is first-order accurate in time. We compare our 
results with the monolithic scheme of Badia, Quaini, and Quarteroni [HI E] 
showing excellent agreement and comparable accuracy. 

In hemodynamics, the coupling between fluid and structure is highly non- 
linear due to the fact that the fluid and structure densities are roughly the 
same, making the inertia of the fluid and structure roughly equal. In this 
regime, classical partitioned loosely coupled (or explicit) numerical schemes, 
which are based on the fluid and structure sub-solvers, have been shown to be 
intrinsically unstable [7J due to the miss-match between the discrete energy 
dictated by the numerical scheme, and the continuous energy of the coupled 
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problem. This has been associated with the explicit role of the "added mass 
effect", introduced and studied in [7J. To rectify this problem, the fluid and 
structure sub-solvers need to be iterated until the energy balance at the dis- 
crete level approximates well the energy of the continuous coupled problem. 
The resulting strongly coupled partitioned scheme, however, gives rise to 
extremely high computational costs. 

To get around these difficulties, several different loosely coupled algo- 
rithms have been proposed that modify the classical strategy in coupling the 
fluid and structure sub-solvers. The method proposed in [10] uses a simple 
membrane model for the structure that can be easily embedded into the fluid 
equations and appears as a generalized Robin boundary condition. In this 
way the original problem reduces to a sequence of fluid problems with a gen- 
eralized Robin boundary condition that can be solved using only the fluid 
solver. A similar approach was proposed in [11], where the fluid and struc- 
ture are split in the classical way, but the fluid and structure sub-problems 
were linked via novel transmission (coupling) conditions that improve the 
convergence rate. Namely, a linear combination of the dynamic and kine- 
matic interface conditions was used to artificially redistribute the fluid stress 
on the interface, thereby avoiding the difficulty associated with the added 
mass effect. 

A different stabilization of the loosely coupled (explicit) schemes was pro- 
posed in [T2] which is based on Nitsche's method [33] with a time penalty term 
giving L 2 -control on the fluid force variations at the interface. We further 
mention the scheme proposed in [T3] , where Robin- Robin type preconditioner 
is combined with Krylov iterations for the solution of the interface system. 

For completeness, we also mention several semi-implicit schemes. The 
schemes proposed in [JSJ [THl HZ] separate the computation of fluid velocity 
from the coupled pressure-structure velocity system, thereby reducing the 
computational costs. Similar schemes, derived from algebraic splitting, were 
proposed in J9j [18]. We also mention [19] where an optimization problem is 
solved at each time-step to achieve continuity of stresses and continuity of 
velocity at the interface. 

In our work we deal with the problems associated with the added mass 
effect by: (1) employing the kinematic coupling condition implicitly in all 
the sub-steps of the splitting, as in the kinematically coupled scheme first 
introduced in ]5]; (2) treating the fluid sub-problem together with the vis- 
cous part of the structure equations so that the structure inertia appears 
in the fluid sub-problem (made possible by the kinematic coupling condi- 
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tion), giving rise to the energy estimates that mimic those in the continuous 
problem. In this step, a portion of the fluid stress and the viscous part of 
the structure equations are coupled weakly, and implicitly, thereby adding 
dissipative effects to the fluid solver and contributing to the overall stability 
of the scheme (although the scheme is stable even if viscoelasticity of the 
structure is neglected). The modification of the Lie splitting introduced in 
this manuscript uses the remaining portion of the normal fluid stress (the 
pressure) to explicitly load the structure in the elastodynamics equations, 
significantly increasing the accuracy of our scheme when compared with the 
classical kinematically coupled scheme |3J, and making it comparable to that 
of the monolithic scheme presented in [HI E] • 




Figure 1: Longitudinal displacement in a carotid artery measured using in vivo ultrasound 
speckle tracking method. The thin red line located at the intimal layer of the arterial wall 
shows the direction and magnitude of the displacement vector, showing equal magnitude 
in longitudinal and radial components of the displacement [20] . 

To deal with the motion of the fluid domain, we implemented an Arbitrary 
Lagrangian-Eulerian (ALE) approach. In addition to the ALE method [2T1 
E21 [231 EH HH1 E5l EHl E3 EHJ , the Immersed Boundary Method [291 EH EH E2 
|33| 134"] has been very popular in problems with moving domains, especially 
when the structure is completely immersed in the fluid domain. We also 
mention the Fictitious Domain Method combined with the mortar element 
method or ALE method [351 [36] , the Lattice Boltzmann method [374 EH EHl 
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ED], the Coupled Momentum Method [H], and the Level Set Method g2]. 

Even though other viscoelastic models have been proposed in literature 
to study FSI between blood and arterial walls, see e.g., [US HH SSI EE] , to the 
best of our knowledge, the present manuscript is the first in which both radial 
and longitudinal displacement of a thin, viscoelastic structure are captured. 
The main motivation for this work comes from recent developments in ul- 
trasound speckle tracking techniques (see [47], Chapter 8), which enabled in 
vivo measurements of both the longitudinal and diameter vessel wall changes 
over the cardiac cycle, indicating that longitudinal wall displacements can be 
comparable to the radial displacements, and should be included when study- 
ing tissue movement [1EJ El EE] • See Figure [I] This is particularly pronounced 
under adrenaline conditions during which the longitudinal displacement of 
the intima-media complex increases by 200%, and becomes twice the magni- 
tude of radial displacement [4~9] . 

The structure model capturing both radial and longitudinal displacement 
is presented next. 

2. The cylindrical linearly viscoelastic Koiter shell model 

In this section we present the main steps in the derivation of the cylin- 
drical linearly viscoelastic Koiter shell model that includes both longitudinal 
and radial components of the displacement. 

Consider a clamped cylindrical shell of thickness h, length L, and refer- 
ence radius of the middle surface equal to R. See Figure [2j This reference 
configuration of the thin cylindrical shell will be denoted by 

T = {x = (R cos 9, R sin 9, z) E R 3 : 9 E (0,2w),z E (0,L)}. (1) 

Displacement of the shell corresponds to the displacement of the shell's 




Figure 2: Left: Cylindrical shell in reference configuration with middle surface radius R 
and shell thickness h. Right: Deformed shell. 
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middle surface. Following the basic assumptions under which the Koiter 
shell model is valid [50], we assume that the walls of the cylinder are thin 
with respect to its radius (i.e., h/R <C 1), homogeneous, and deform linearly 
(i.e., the displacement and the displacement gradient are both small). We 
will be assuming that the load exerted onto the shell is axially symmetric, 
leading to the axially symmetric displacements, so that the displacement in 
the direction is zero, and nothing in the problem depends on 9. Thus, 
displacement r? will have two components, the axial component r] z and the 
radial component r\ r . 

To introduce the viscous effects into the linearly elastic Koiter shell model, 
we must assume that displacement is a function of both position and time. 
Thus ri(z,t) = (r] z (z,t),r] r (z,t)). The derivatives with respect to the spatial 
variable will be denoted by r?', and with respect to the temporal variable by 
77. 

To define the Koiter shell model we need to define the geometry of de- 
formation, and the physics of the problem which will be described by the 
dynamic equilibrium equations (the Newton's second law of motion). 

Geometry of Deformation. The axially symmetric configuration of a 
Koiter shell can suffer stretching of the middle surface, and flexure (bending). 
The stretching of the middle surface is measured by the change of metric 
tensor, while flexure is measured by the change of curvature tensor. The 
change of metric and the change of curvature tensors for a cylindrical shell 
are given, respectively by [51] 



71*7 J 















Rr} r 


Tj r _ 



(2) 



The dynamic equilibrium equations. The total energy of the linearly 
elastic Koiter shell is given by the sum of contributions due to stretching 
and flexure. The corresponding weak formulation will thus account for the 
internal (stretching) force and bending moment. For a linearly viscoelastic 
Koiter shell each of the contributions will also have an additional component 
to account for a viscoelastic response of the Koiter shell. The elastic response 
will be modeled via the elasticity tensor A, while the viscous response will 
be modeled by the viscosity tensor B. Viscoelasticity will be modeled by the 
Kelvin- Voigt model in which the total stress is linearly proportional to strain 
and to the time-derivative of strain. 
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More precisely, introduce the elasticity tensor A, [51], to be denned as: 



AE 



where A r 



2Ea 
I -a 2 

1 



f A c • E)A C 



2E 
1 + cr 



A C EA C , E G Sym(I 



(3) 



is the first fundamental form of the middle surface, 



A c 



1 



R 2 

is its inverse, and ■ denotes the scalar product 
A - B :=Tr(AB T ), A,BeM 2 (R). 



Here E is the Young's modulus and a is the Poisson's ratio. 
Similarly, define the viscosity tensor B by: 



BE 



2E v a v 



A c • E)A C 



2R 



1 + cr. 



— A C EA C , E G Sym(I 



(4) 



Here E v and cr„ correspond to the viscous counterparts of the Young's mod- 
ulus E and Poisson's ratio a. Then, for a linearly viscoelastic Koiter shell 
model we define the internal (stretching) force 



and bending moment 



N := ^ 7 (T7) + ^ 7 (t)), 



M := ^(77) + ^Bg(fj). 



(5) 



(6) 



The weak formulation of the linearly viscoelastic Koiter shell is then given 
by the following: for each t > find r](-,t) G V c such that V£ G V c 



k f {A 1 {r 1 ) + B 1 {r } ))- 1 {£)Rdz+ 1 ^ 



(AQ(<n)+BQ(ri))-Q(£)Rdz 



+p s h 



d r] 



£Rdz 



f ■ £Rdz, 



where p s denotes the volume shell density and 
V c = H*(0,L)xH*(0,L) 



(7) 



{£ = (Z z ,Zr)eH\0,L) : ^)^(L) = 0,e:(0)=e(L)=0}. 
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The first term on the left-hand side of ([7]) multiplying h/2 captures the 
membrane effects, while the second term on the left hand-side multiplying 
h 3 /24 captures the flexural effects of the Koiter shell. 

Components of the forcing term f = (f z , f r ) T are the surface densities in 
the reference configuration of the axial and radial force. The corresponding 
dynamic equilibrium equations in differential form can be written as follows: 

d 2 7] r d 2 7] r 8t] z d^T\ r drj r d 3 7] r 



2 



dt 2 ' ' dz 2 dz dz 4 dt dtdz 



where 



hE h 2 h 3 Ea h Ea 

^0 — -fvTT-x 9\K l + TTTTv? J; ^1 — —-^rrz ?rr> ^: 



R 2 {l-a 2 y 12B 2 " 1 6R 2 (l-a 2 )' z Rl-a 2 
_ hE h 3 E 

C3 "T^ 2 "' U "121^2' 

h 



(10) 



and 



t-j E v <j v 



<Jz 1 — or 



The boundary conditions for a clamped Koiter shell problem are given by 

r}(0,t) = ri(L,t) = 0, ^(0,0 = ^(^,0 = 0. (11) 
For a mathematical justification of the Koiter shell model please see [HJJ E21 



3. The fluid-structure interaction problem 

We consider the flow of an incompressible, viscous fluid in a two-dimensional 
channel of reference length L, and reference width 2R, see Figure [3} The lat- 
eral boundary of the channel is bounded by a thin, deformable wall, modeled 
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by the linearly viscoelastic Koiter shell model, described in the previous sec- 
tion. We are interested in simulating a pressure-driven flow through the 
deformable 2D channel with a two-way coupling between the fluid and struc- 
ture. Without loss of generality, we consider only the upper half of the fluid 
domain supplemented by a symmetry condition at the axis of symmetry. 
Thus, the reference domain in our problem is given by 

Q Q ■= {( z , r)|0 < z < L, < r < R}. 

Here z and r denote the horizontal and vertical Cartesian coordinates, re- 
spectively. See Figure [3j 

Remark 1. It is worth mentioning here that while the fluid flow will be 
modeled in 2D, the thin structure equations, described in the previous sec- 
tion, are given in terms of cylindrical coordinates, assuming axial symmetry. 
It is standard practice in 2D fluid-structure interaction studies to use thin 
structure equations that are derived assuming cylindrical geometry. This is 
because cylindrical structure models account for the circumferential stress 
that "keeps" the top and bottom boundary of the structure "coupled to- 
gether" when they are loaded by the stresses exerted by the fluid, thereby 
giving rise to physiologically reasonable solutions. 

The mathematical model for the corresponding fluid-structure interaction 
problem can be defined as follows. The fluid domain, which depends on 
time, is not known a priori. The location of the lateral boundary, defined in 
Lagrangian framework, is given by T(t) = {(z + r] z (z,t),R + r] r (z,t)) \ z G 
(0,L)} for t G (0,T). Throughout the rest of the manuscript we will be 
denoting the Lagrangian coordinates by x = (z,f). The displacement of the 
boundary will always be given in Lagrangian framework. However, we will 
omit the hat notation on 17 for simplicity. 

We will be assuming that for each t G (0, T) the boundary of the fluid 
domain is Lipschitz continuous (see Figure [3]), and that it's lateral boundary, 
in Eulerian framework, can be described by a Lipschitz continuous function 

g(- ; t) : (0, L) — >■ E, g{- ;t) : z H> g(z;t), for each t G (0,T) 

so that, in Eulerian framework, 

T(t) = {(z,g(z;t)),ze(0,L)}, for te(0,T). 



9 



The fluid domain is given by 

Q(t) = {{z,r) G R 2 ; < z < L, < r < g{z;t)}, for ie(0,T). (12) 

The inlet boundary will be denoted by r in , the outlet boundary by r out , the 
symmetry (bottom) boundary for which r = by T , so that 

0ft(t) = r in u r(t) u r out u r . 





r 









R 

Tin 












To 


L 2 













Figure 3: Deformed domain Cl(t). 



The flow of a viscous, incompressible, Newtonian fluid is governed by the 
Navier-Stokes equations 

P/^ + u.Vuj = V-<r in Cl(t) for t G (0, T), (13) 
V-u = in fi(£) for t G (0, T), (14) 

where u = (u z ,u r ) is the fluid velocity, p is the fluid pressure, pf is the fluid 
density, and a is the fluid stress tensor. For a Newtonian fluid the stress 
tensor is given by cr = —pi + 2/iD(u), where /i is the fluid viscosity and 
D(u) = (Vu + (Vu) r )/2 is the rate-of-strain tensor. 

At the inlet and outlet boundary we prescribe the normal stress: 

<rn| in (0,r,t) = -p in (t)n\ in on (0, R) x (0, T), (15) 
an\ out (L,r,t) = -p out {t)n\ out on (0, R) x (0, T), (16) 

where n in /ri out are the outward normals to the inlet/outlet boundaries, re- 
spectively. These boundary conditions are common in blood flow model- 
ing uniESlES]. 

At the bottom boundary r = we impose the symmetry conditions: 

-t^(z, 0, t) = 0, u r (z, 0, t) = on (0, L) x (0, T). (17) 



10 



The upper boundary T(t) represents the deformable channel wall, whose 
dynamics is modeled by ([8])-(|9|. The structure equations are supplemented 
with the following boundary conditions 

r)(0,t) — r](L,t) — d z r) r (0,t) — d z r) r (L,t) — on(0,T). (18) 

Initially, the fluid and the structure are assumed to be at rest, with zero 
displacement from the reference configuration 

u = 0, 77 = 0, ^ = 0. (19) 

The fluid and structure are coupled via the kinematic and dynamic bound- 
ary conditions [57]: 

• Kinematic coupling condition describes continuity of velocity 

u(z + Vz (z,t),R + ri r (z,t),t) = ^(z,t) on (0, L) x (0,T). (20) 

• Dynamic coupling condition describes balance of contact forces, 
namely, it says that the contact force exerted by the fluid is equal but 
of opposite sign to the contact force exerted by the structure to the 
fluid: 

U = -J ah\ m -e,on{0,L) x (0,T), (21) 
f r = -J ah\ m -e r on(0,L) x (0,T), (22) 

where 




J=\ i + ^ + PSr ( 23 ) 



denotes the Jacobian of transformation from the Eulerian to the La- 
grangian framework, and ern denotes the normal fluid stress on the 
reference domain Q = (0,L) x (0,-R). Here e z = (1,0) and e r = (0, 1) 
are the standard unit basis vectors, and n is the outward normal to 
the deformed domain. 
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3.1. The energy of the coupled F SI problem 

To formally derive the energy of the coupled FSI problem we multiply 
the structure equations by the structure velocity, the balance of momen- 
tum in the fluid equations by the fluid velocity, integrate by parts over the 
respective domains using the incompressibility condition, and add the two 
equations together. The dynamic and kinematic coupling conditions are then 
used to couple the fluid and structure sub-problems. The resulting equation 
represents the total energy of the problem. 

We start by first considering the Koiter shell model for the structure. 
We recall the weak formulation of the clamped Koiter shell given by ([7]). In 
([7]) we replace the test function £ by the structure velocity ^ and integrate 
by parts over (0, L) to obtain the following energy equality of the clamped 
Koiter shell: 



+ '-. 



l+cr 



+ 



III 
24 



l+cr 



E v 



l+cr v 



d_ J p s h 
dt] 2 



dr\ z 



+ 



L 2 (0,L) 



dt 



1+CT 



+ 



Psh 



L 2 (0,L) 



drjr 
dt 



L 2 (0,L) 



djQz 



dz 



+ 



Ecr 



L 2 (0,L) 



1-cr 2 



dr) z _|_ rjr 
dz ' R 



+ 



L 2 (0,L) 
2 



1 + CT 



d 2 r] r 



drjr 



Rdt 



+ 



Ey 



L 2 (0,L) 



l + cr„ 



d 2 y z 



L 2 (0,L) 
2 



I _Ecr 

h 1-a 2 



L 2 (0,L) 
2 



dzdt 



L 2 (0,L) 



E v a v 



+ ^2 



d 2 rj r _|_ rjr_ 
dz 2 ^ R 2 



d 2 r) z 1 dr/r 
dzdt ' Rdt 



L 2 (0,L) 
2 



+ 



hi 
21 



E v 



l+cr„ 



R 2 dt 



+ 



E v 



L 2 (0,L) 



1 + 01; 



9 3 



')r 



dz 2 dt 



_ E v cr v 
h 1-oi 



L 2 (0,L) 

Jo f ' ~bt dz- 



L 2 (0,L) 
2 



d 3 r] r 1 dr) r 
dz 2 dt R2dt 



L 2 (0,L) 



(24) 



The first two terms under the time-derivative correspond to the kinetic energy 
of the Koiter shell. The terms in the second and third row correspond to the 
elastic energy of the Koiter shell (the terms multiplying h are the membrane 
energy, while the terms multiplying h 3 correspond to the flexural (bending) 
energy). The terms in the fourth and fifth row correspond to the viscous 
energy of the viscoelastic Koiter shell, while the last term corresponds to the 
work done by the external loading, which comes form the fluid stress via the 
dynamic coupling condition. 
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To deal with the fluid sub-problem we multiply the momentum equation 
in the Navier-Stokes equations by u and integrate by parts over Q(t), using 
the incompressibility condition along the way. With the help of the following 
identities 



du 

~dt 



udx 



{u ■ V)u ■ udx 



'n(t) 
one obtains 



1 d 

1 

2 



| IX | dx — ; 

it| 2 it ■ ndS 1 , 



I it I u ■ ndS 1 , 



1 d 

2dt 



pf\\ u \\mm) \ + M\ d (. u )\\l* 



(«(*)) 



Pin(*)w*|*=o^+ / Pout(t)u z \ z=L dr = / crn-UGLS" 
o </o Jr(t) 

The integral on the right-hand side can be written in Lagrangian coordinates 
as 

/ trn-ud3 = [em u] \(i +v 4z,t),R+vr(B,t)) J dz (25) 
</r(t) Jo 

where J is the Jacobian of transformation from the Eulerian to Lagrangian 



framework, given by (|23|). Now we use the kinematic and dynamic lateral 

' L drj 



boundary conditions (20)-(22) to obtain 

-L 



I [crn ■ u) \(z+- nz (z,t),R+ Vr (z,t)) J dz — — f 
Jo Jo 



at 



dz. 



(26) 



Thus, the fluid sub-problem coupled with the structure satisfies 



1 d 
2 



I P/IHIW)) \+2»\\D(u)\\l 2 



R pR 

Pin(t)u z \ z=0 dr + / Pout(t)u z \ z=L dr 



L 



f 9l l A" 



(27) 



By adding (24) and (27), the right-hand sides of the two equations cancel 
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out and one obtains the energy equality for the FSI problem: 
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Fluid Viscous Energy 



(28) 



The coefficients E v ,a v are defined after equation Q. Therefore, we have 
shown that if a solution to the coupled fluid-structure interaction problem (fsh 



- (22) exists, then it satisfies the energy equality (28). This equality says that 
the rate of change of the kinetic energy of the fluid, the kinetic energy of the 
structure, and the elastic energy of the structure, plus the viscous energy of 
the structure, plus the viscous energy of the fluid, is equal to the work done 
by the inlet and outlet data. 



4. The numerical scheme 



To solve the fluid-structure interaction problem (|8j)-(22), we propose an 
extension of a loosely coupled partitioned scheme, called the kinematically 
coupled scheme, first introduced in [5]. 

The classical kinematically coupled scheme introduced in [5] is based on 
a time-splitting approach known as the Lie splitting [58]. The viscoelas- 
tic structure is split into its elastic part and the viscous part. The viscous 
(parabolic) part is treated together with the fluid, while the elastic (hyper- 
bolic) part is treated separately. The inclusion of the viscous part of the 
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structure into the fluid solver as a boundary condition in the weak formula- 
tion accounts for the fluid load onto the structure as it simultaneously deals 
with the "added mass effect" [7j. This adds dissipative effects to the fluid 
solver that contribute to the stability of the scheme. This approach provides 
a desirable discrete energy inequality making this scheme stable even when 
the density of the fluid is equal to the density of the structure, which is the 
case in the blood flow application. The elastic part of the structure, which is 
solved separately, communicates with the fluid only via the kinematic cou- 
pling condition in the classical kinematically coupled scheme. The fluid stress 
does not appear in this step, as it is used as a loading to the viscous part of 
the structure in the weak formulation of the fluid sub-problem. 

In this manuscript we will change this approach by additionally splitting 
the normal stress into a fraction that loads the viscous part of the structure, 
and a fraction (pressure) that loads the elastic part of the structure. This 
splitting is done using a modification of the Lie splitting scheme in a way 
which significantly increasases accuracy. 

Thus, in this manuscript, the kinematically coupled scheme is extended 
and improved to achieve the following two goals: 

1. Capture both the radial and longitudinal displacement of the linearly 
viscoelastic Koiter shell for the underlying fluid-structure interaction 
problem. 

2. Increase the accuracy of the kinematically coupled scheme by introduc- 
ing a new splitting strategy based on a modified Lie's scheme. 

This version of the kinematically coupled scheme retains all the advan- 
tages of the original scheme, which include: 

• The scheme does not require sub-iterations between the fluid and struc- 
ture sub-solvers to achieve stability. 

• The scheme is modular, allowing the use of one's favorite fluid or struc- 
ture solvers independently. The solvers communicate through the ini- 
tial conditions. 

• Except for the pressure, the fluid stress at the boundary does not need 
to be calculated explicitly. 
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4-1. The Lie scheme 

To apply the Lie splitting scheme the problem must first be written as a 
first-order system in time: 

d ^ + A{<P) = 0, in(0,T), (29) 
0(0) = O , (30) 

where A is an operator from a Hilbert space into itself. Operator A is then 
split, in a non-trivial decomposition, as 



A = J2 A >- (31) 



i=l 



The Lie scheme consists of the following. Let At > be a time discretization 
step. Denote t n = nAt and let <p n be an approximation of <f)(t n ). Set <f)° = <po- 
Then, for n > compute <p n+1 by solving 

t^i + Aifc) = in(f\r +1 ), (32) 
(f) t {t n ) = n+ ( i - 1 )/ / , (33) 

and then set (j) n+i/I = 4 (t n+1 ), for % = 1, . . . .1. 



This method is first-order accurate in time. More precisely, if (29) is 
defined on a finite-dimensional space, and if the operators Ai are smooth 
enough, then ||</>(t n ) — n || = O(At) |58]- In our case, operator A that is 



associated with problem (|8j)-(22) will be split into a sum of three operators: 

1. The time-dependent Stokes problem with suitable boundary conditions 
involving structure velocity and fluid stress at the boundary. 

2. The fluid advection problem. 

3. The elastodynamics problem for the structure loaded by the fluid pres- 
sure. 

These sub-problems are coupled via the kinematic coupling condition and 
via fluid pressure appearing in the elastodynamics problem. The kinematic 



coupling condition also plays a key role in writing problem rt8h-( 22 ) as a 



first-order system, based on which the Lie splitting can be performed. 
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4-2. The first-order system in ALE framework 

To deal with the motion of the fluid domain we adopt the Arbitrary 
Lagrangian-Eulerian (ALE) approach j2TJ E2J [56]. In the context of finite 
element method approximation of moving-boundary problems, ALE method 
deals efficiently with the deformation of the mesh, especially near the in- 
terface between the fluid and the structure, and with the issues related to 
the approximation of the time-derivatives du/dt ~ (u(t n+1 ) — u(t n ))/At 
which, due to the fact that Cl(t) depends on time, is not well defined since 
the values u(t n+1 ) and u(t n ) correspond to the values of u defined at two 
different domains. ALE approach is based on introducing a family of (arbi- 
trary, invertible, smooth) mappings At defined on a single, fixed, reference 
domain Cl such that, for each t £ (t ,T), At maps the reference domain 
Cl = (0,L) x (0,-R) into the current domain Cl(t) (see Figure 4J): 



At : Cl C M 2 -»■ Cl(t) C 



x = A(x) £ Cl(t), for x £ Cl. 



In our approach, we define At to be a harmonic extension of the mapping 




Figure 4: At maps the reference domain fl into the current domain f2(t). 



g that maps the boundary of Cl to the boundary of Cl(t) for a given time 
t. More precisely, in our case Cl := (0, L) x (0,-R), and so At is a harmonic 
extension of 

g : dCl ->■ dCl(t) 
onto the whole domain Cl, for a given t : 



At 



AA t 

\dCl\t 



in O, 

9- 
0. 



To rewrite system l8[)-( 22 ) in the ALE framework we notice that for a 



function / = /(x, t) defined on fi(t) x (0, T) the corresponding function 
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/ := / o At defined onfix (0, T) is given by 

f(x,t) = f(A t (x),t). 
Differentiation with respect to time, after using the chain rule, gives 

Of 



df 
dt 



dt 



+ w-V/, 



(34) 



where w denotes domain velocity given by 

dAAx) 



w(x,t) 



dt 



(35) 



To write syste m (|9[ )-(22) in first-order form, we utilize the kinematic 
coupling condition ( |20[ ). Written in the ALE framework, our problem now 
reads: Find u = (u z ,u r ), rj = (r] z ,r] r ), with u(x, t) = u(^4 t (x),i) and u\f = 
u(z, R, t), such that 



Pf 



+(u-w)-Vuj = V-o", inQ(t) x (0,T), (36) 
V-u = infl(t) x (0,T), (37) 
with the kinematic and dynamic coupling conditions holding on T(t): 
drj 



dt 
p s h 



u 



on (0,L) x (0,T), 
dr] r d 2 rj z 



(38) 



dt 



^2^. C3 o n 



<9(u 



(9,2 



zirv 



1 + 



drjg 
dz 



+ 



drjr 
dz 



ah\ r{ty e z on (0,L) x (0,T),(39) 



Psh 



-Dr 



dt 

d 2 (u 



r\v> . ^ ^ «9V dr) z d\ A 



r Ipv 



fe 2 



<9(w 



zirv 



+ £>4 



<9 4 fw 



<9z 4 



1 + 



<9r^ 
dz 



+ o^lrft) ■ e r on (0, L) x (0, T), (40) 
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and the following boundary conditions on r in U r out U r : 

du 7 



Or 



(z,0,t) = u r (z,0,t) = on r , 



u(0,R,t) = u(L,R,t) = 0, f]\ z=0 ,L 



drj r 
dz 



=0,L 



crn| in (0,r,t) = -j?» n (*)n| in , 

crn| out (L,r,t) = -p ou t(t)n| out on (0,i?) x (0,T). 

At time t — the following initial conditions are prescribed: 



u| t=0 = 0, 77| t=0 = 0, 



drj 



0. 



(41) 

(42) 

(43) 
(44) 

(45) 



t=o 



Notice how the kinematic coupling condition is used in (39) and (40) to 
rewrite the viscous part of the structure equations in terms of the trace 
of the fluid velocity on T(t). This will be used in the splitting algorithm 
described below. 



Remark 2. As shown in [59], if we discretise (35) as 

Ar(x) 



x 



W[X,T 



+ 0(t) 



we obtain a linear affine transformation for A T 

Ar(x) — X + TW(X, t) + 0(t). 



(46) 



(47) 



It can be easily shown that, using this transformation, spatial partial deriva- 
tives of a function on a domain O(r) are equal to the derivatives of the same 
function on the reference domain Cl, plus an error 0(t) [59]. We avoid deal- 
ing with this problem by writing only the time- derivative on the reference 
domain, and leaving the spatial derivatives evaluated on the current domain. 

4-3. Details of the operator- splitting scheme 

We split the first-order system (36)-(45) into three sub-problems. The 
fluid problem will be split into its viscous part and the pure advection part 
(incorporating the fluid and ALE advection simultaneously). The fluid stress 
efn will be split into two parts, Part I and Part II: 



(fn = an + (3pn — /3pn, 
(i) in) 



(48) 
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where (3 is a number between and 1, < /3 < 1, with f3 = corresponding 
to the splitting introduced in [5]. As will be shown later, the accuracy of the 
scheme changes as the value of /3 increases from to 1. The numerical results 
presented in this manuscript will correspond to the value of /3 = 1, since our 
numerical investigation showed that /3 = 1 provides highest accuracy for 
Examples 1 and 2 presented in this manuscript. 

The viscoelastic structure equations will be split into their viscous part 
and the elastic part. These are combined into a splitting algorithm in the 
following three steps. 

• Step 1. Step 1 involves solving the time-dependent Stokes prob- 
lem, incorporating the viscous part of the structure and Part I of the 
fluid stress via a Robin-type boundary condition. The time-dependent 
Stokes problem is solved on a fixed domain Q(t n ). The problem reads 
as follows: 

Find u,p and rj, with u(x, t) = u(^4 t (x),t) such that for t G (t n ,t n+1 ), 
with p n and r) n obtained at the previous time step: 

' PfW L = V<T ' V-u = infi(r) 

fOM) = on(0,L) x (f\ 

Ps h 9 -^ - D 2 9 -^ - D 3 %M + V( 1 + ¥) 2 + (¥) 2 (^)k") • e z 
= -\/(l + ¥) 2 + (¥) 2 (^)lr(^) • e z on (0, L) x (r,^ 1 ), 

Psh 9 -^ + D u r \ t - D t q^ + D 2 9 -^ + D/-^ 
+/ 3 V / ( 1 + ¥) 2 + (¥) 2 (^)lr( i ") " e r 
= -V(l + ¥) 2 + (¥) 2 (^)lr(^) • e r on (0, L) x 

with the following boundary conditions on Ti n U r out U To: 
— {z,0,t)= u r {z,0,t) = onr , 

u(0, R,t) = u{L,R,t) = 0, 
crn\ in = -p in (t)n\ in on r in , crn\ out = -p ou t{t)n\ out on r out , 
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and initial conditions 



u(t n ) = u n , n(t n ) = rf. 

Then set u n+1 / 3 = u(r +1 ), r] n+1 ' z = r]{t n+1 ), p n+1 = p(t n+1 ). 
Note that here we used only Part I of the fluid stress. 

• Step 2: Solve the fluid and ALE advection sub-problem defined on a 
fixed domain Q(t n ). The problem reads: Find u and 77 with u(x, t) = 
u(A(x),t), such that for t E (t n ,t n+1 ) 

' §f | 4 + (u n+1 / 3 - w n+1 / 3 ) • Vu = 0, in tt(t n ) 

< §fOM) = o on(o,L)x(r,r +1 ), 
k M-^ = o, on (o,l) x (r,r +1 ), 

with boundary conditions: 

u = u n+1 / 3 on T^ +1/3 , where 

pn+l/3 = { x G R 2| x G (u n+l/3 _ w n+l/3) . n < Q}, 

and initial conditions 

u(r) = u n+1/3 , »7(i n ) = r7 n+1/3 . 
Then set u n+2 / 3 = u(r +1 ), r] n+2 ^ = r]{t n+l ). 

• Step 3: Step 3 involves solving the elastodynamics problem for the 
location of the deformable boundary by involving the elastic part of 
the structure which is loaded by Part II of the normal fluid stress. 
Additionally, the fluid and structure communicate via the kinematic 
lateral boundary condition which gives the velocity of the structure in 
terms of the trace of the fluid velocity, taken initially to be the value 
from the previous step. The problem reads: Find u and 77, with p n+1 
computed in Step 1 and rf 1 obtained at the previous time step, such 
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that for t E (t n ,t n+1 ) 
= 0, in Q(r) 

%{z,t)=u\ t on(0,L)x(r,r +1 ), 



( du 

at 



^ = ^ 1 + S) 2 + (¥) 2 (^ n )lr(*") • e z on (0, L) x (r,^ 1 ), 

Ps n — t r O ^ r - ^1-^5" + ^2g7 + Wg^r 

= /37( 1 + ^) 2 + (^) 2 (P^ n )lr(t-) ' e r on (0,L) x 



with boundary conditions: 

I _ dr ^n _ n 

»7U=o,l — -q^\z=o,l — u; 

and initial conditions: 

u(t n ) = u n+2 / 3 , f7(i n ) = r7 n+2 / 3 . 

Then set u n+1 = u(r +1 ), r] n+1 = t]{t n+l ). 
Do t n = t n+1 and return to Step 1. 

Remark 3. Note that the outward normal to the lateral boundary can be 
written as 

Using this equality, we can take n in Step 3 implicitly, which upon substi- 
tuting u|p by ^ leads to the following system 

- (c 2 - - = o on (o,l) x 

+ Coifr - + (c 2 - m^-t + 

= ffi$ t) on(0,L)x(r,r +1 ), 
where p n+1 is pressure computed in Step 1. 

Remark 4. The trace of the pressure, used in Step 3 to load the structure, 
needs to be well-defined. In general, one expects the pressure for a Dirichlet 
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problem defined on a Lipschitz domain to be in L 2 (Q), which is not sufficient. 
Several works, see e.g., [60.[ EH E2] , indicate that, under certain compatibility 
conditions, the solution of the related class of moving boundary problems has 
higher regularity, allowing the definition of the trace of the pressure on the 
moving boundary. In fact, we can show that for our problem, under certain 
compatibility conditions at the corners of the domain, the pressure belongs 
to W lfi / 7 {n), which is more than sufficient for the trace to be well-defined 
on r(r) [63]. 



5. Numerical results 

We present three examples that show the behavior of our scheme for dif- 
ferent parameter values. Example 1 below, corresponds to the benchmark 
problem suggested by Formaggia et al. in [M] to study the behavior of FSI 
scheme for blood flow. The structure model in this case is a linearly viscoelas- 
tic string model capturing only radial displacement, with the coefficient that 
describes bending rigidity given in terms of the shear modulus and Timo- 
shenko correction factor, which is different from the corresponding Koiter 
shell model coefficient. The structural viscosity constant and the Youngs 
modulus in this model are both relatively small. 

In Example lb we supplemented this model by an equation describing dy- 
namics of longitudinal displacement, obtained from the Koiter shell model. 
The corresponding equation for longitudinal displacement is degenerate, as it 
does not involve any spatial derivatives as in the equation for radial displace- 
ment, and there are no coupling terms between the radial and longitudinal 
displacement. 

In Example 2 we consider the full Koiter shell model capturing both ra- 
dial and longitudinal displacement, and the coupling between the two. The 
coefficients in the model are given by those associated with the derivation of 



the Koiter shell, see (10). The values of Youngs modulus, Poisson ratio, and 
shell thickness are the same as in Examples 1 and lb, with the structural 
viscosity coefficients small, and related to the one in Example 1. This exam- 
ple could be used as a benchmark problem for testing numerical methods for 
FSI in which the structure is modeled as a linearly viscoelastic Koiter shell. 

The last example concerns simulation of blood flow through a section of 
the Common Carotid Artery (CAA). Numerical simulations were compared 
with experimental data showing very good agreement. 
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The following values of fluid and structure parameters are used in Exam- 
ples 1 and 2, listed below: 



Parameters 


Values 


Parameters 


Values 


Radius R (cm) 


0.5 


Length L (cm) 


6 


Fluid density pj (g/cm 3 ) 


1 


Dyn. viscosity p (poise) 


0.035 


Wall density p s (g/cm 3 ) 


1.1 


Wall thickness h s (cm) 


0.1 


Young's mod. E (dynes/cm 2 ) 


0.75 x 10 e 


Poisson's ratio a 


0.5 



Table 1: Geometric parameters, and fluid and structural parameters that are used in 
Examples 1 and 2 presented in this section. 



Parameter (3, introduced in (48), which appears in Step 1 and Step 3 
of our numerical scheme, can vary between and 1, where the value of 
= corresponds to the kinematically coupled scheme presented in [5] . The 
change in j3 is associated with the change in the accuracy of the scheme (not 
the stability). For the test cases presented in Examples 1 and 2, the value 
of (3 = 1 provides the highest accuracy. We believe that the main reason 
for the gain in accuracy at (3 = 1 is the strong coupling between the fluid 
pressure (which incorporates the leading effect of the fluid loading onto the 
structure) and structure elastodynamics, which is established for (3 = 1 in 
Step 3 of the splitting, described above. Namely, in the case (3 = 1, the 
structure elastodynamics problem is forced by the entire contribution of the 
fluid pressure, which appears in Step 3 as a loading onto the sturucture. The 
elastodynamics of the structural problem is, therefore, directly coupled to 
the entire fluid pressure for (3 = 1. It remains to be investigated how does 
the optimum choice of (3 (providing highest accuracy for a given problem), 
depend on the coefficients in the problem. 

5.1. Example 1: The benchmark problem with only radial displacement. 

We consider the classical FSI test problem proposed by Formaggia et 
al. in [M]. This problem has been used in several works as a benchmark 
problem for testing the results of fluid-structure interaction algorithms for 
blood flow [TQ1 EEl El El E] • The structure model for this benchmark problem 
is of the form 

Psh^ - kGh— + 7^ = /, (50) 
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with absorbing boundary conditions at the inlet and outlet boundaries: 



at z = (51) 



at z = L. (52) 
Here G = on E . , is the shear modulus and k is the Timoshenko shear correc- 

2(1+17) 

tion factor. The flow is driven by the time-dependent pressure data: 
p in (t) = { 2 L lwJJ ; ^ , PoutW = W G (0,T), (53) 

where p max = 2 x 10 4 (dynes/cm 2 ) and t max = 0.005 (s). The graph of the 
inlet pressure data versus time, is shown in Figure [5] The values of all the 

x 10 4 



T3 

" 

°~ 0.002 0.004 0.006 0.008 0.01 0.012 

time [s] 

Figure 5: The inlet pressure pulse for Examples 1 and 2. The outlet pressure is kept at 0. 
parameters in this model are given in Tables [T] and [2j 



Parameters 


Values 


Shear mod. G(dynes/cm 2 ) 


0.25 x 10 6 


Timoshenko factor k 


1 


Structural viscosity 7 (poise cm) 


0.01 



Table 2: Example 1: Structural parameters considered in Example 1, in addition to those 
listed in Tabic [TJ 

The structural viscosity and Youngs modulus are both very small. For the 
typical physiological values of these parameters see, e.g., [IB]. This means 
that the arterial wall in this example is rather elastic. The relatively large 



dt 



drjr 
dt 



+ 



' kG dr\ r 
Ps dz 



1 kG dr] r 
'., dz 
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value of the coefficient in front of the second-order derivative with respect to z 
(describing bending rigidity), minimizes the oscillations that would normally 
appear in such structures. See Example 2 for more details. 

We implemented this problem in our solver. Since this model captures 
only radial displacement the solver was modified accordingly. The model 
equation (50) can be recovered from the Koiter shell model ^ by setting 



the following values for the coefficients: 

Co — ^2(f~2j) C\ = kGh, C*2 =0, C3 = 0, C4 = 0, 
D =0, D 1 =7, D 2 =0, D 3 =0, D 4 =0. 

The numerical values of these constants are given in Table [3j Homogeneous 



Co 


= 4 x 10 5 


Ci 


= 2.5 x 10 4 


c 2 


= 


c 3 


= 


D 


= 


Dt 


= 10- 2 


D 2 


= 


D 3 


= 



Table 3: Koter shell model coefficients for Example 1. 
Dirichlet boundary conditions for the structure in Step 3 were replaced with 



absorbing boundary conditions (51 )-(52 ). The problem was solved over the 



time interval [0, 0.0 12] s. Propagation of the inlet pressure pulse in terms of 
velocity, displacement, and pressure, vs. time, calculated at the mid-point 
of the tube, are shown in Figure [91 A 2D cartoon of the pressure pulse 



propagating in the tube, is shown in Figure [10 

The numerical results obtained using our modification of the kinemati- 
cally (loosely) coupled scheme proposed in this manuscript, were compared 
with the numerical results obtained using the classical kinematically coupled 
scheme proposed in [3], and the monolithic scheme proposed in [S|. Fig- 
ures [6j [7] and [8] show the comparison between tube diameter, flowrate and 
mean pressure, respectively, at six different times. 

These results were obtained on the same mesh as the one used for the 
monolithic scheme in [8], containing 31 x 11 Fi fluid nodes. More preciesely, 
we used an isoparametric version of the Bercovier-Pironneau element spaces, 
also known as Pi-iso-P2 approximation in which a coarse mesh is used for the 
pressure (mesh size h p ) and a fine mesh for velocity (mesh step h v = hp/2). 

The time step used was At = 10~~ 4 which is the same as the time step 
used for the monolithic scheme, and the kinematically coupled scheme. Due 
to the splitting error, it is well-known that classical splitting schemes usu- 
ally require smaller time step to achieve accuracy comparable to monolithic 
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Figure 6: Example 1: Diameter of the tube computed with the kinematically coupled 
scheme with time step At = 10~ 4 (dash-dot line), implicit scheme used by Quaini in [5] 
with the time step At = 10~ 4 (dashed line) and our scheme with the time step At = 10~ 4 
(solid line). 



schemes. However, the new splitting proposed in this manuscript allows us 
to use the same time step as in the monolithic method, obtaining comparable 
accuracy, as it will be shown next. This is exciting since we obtain the same 
accuracy while retaining the main benefits of the partitioned schemes, such 
as modularity, simple implementation, and low computational costs. 

The kinematically coupled scheme was shown numerically to be first-order 
accurate in time and second order accurate in space [5]. Second-order accu- 
racy in space is retained in the current kinematically coupled /3-scheme, since 
the same spatial discretization of the underlying operators A iy i = 1,2, 3, was 
used in the present manuscript as in |5j. However, due to the new time- 
splitting, the accuracy in time has changed. Indeed, here we show that this 
is the case by studying time-convergence of our scheme. Figure 11 shows a 
comparison between the time convergence of our scheme, the kinematically 
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Figure 7: Example 1: Flowrate computed with the kinematically coupled scheme with 
time step At = 10~ 4 (dash-dot line), implicit scheme used by Quaini in [8 with the time 
step At = 10~ 4 (dashed line) and our scheme with the time step At = 10~ 4 (solid line). 



At 


IIP -Pre/ II 1,2 


L z order 


\\ u - u ref\\ L 2 


L A order 


ll»J - *7re/ILa 


L 2, order 


10" 4 


4.01c + 03 




5.97 




0.003 






(5.65e + 04) 




(136.32) 




(0.0446) 




5 x io _b 


1.57e + 03 


1.35 


4.05 


0.56 


0.0014 


1.1 




(3.36e + 04) 


(0.75) 


(77.91) 


(0.80) 


(0.0264) 


(0.75) 


icr 5 


296.36 


1.04 


1.0 


0.87 


3.17e- 04 


0.92 




(7.27e + 03) 


(0.95) 


(16.27) 


(0.97) 


(0.00576) 


(0.95) 


5 x 10~ b 


134.33 


1.14 


0.46 


1.12 


1.45e-04 


1.13 




(3.3e + 03) 


(1.14) 


(7.36) 


(1.14) 


(0.0026) 


(1.14) 



Table 4: Example 1: Convergence in time calculated at t = 10 ms. The numbers in 
the parenthesis show the convergence rate for the kinematically coupled scheme presented 

coupled scheme, and the monolithic scheme used in [S]. The reference solu- 
tion was defined to be the one obtained with At = 10~ 6 . We calculated the 
absolute L 2 error for the velocity, pressure and displacement between the ref- 
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Figure 8: Example 1: Mean pressure computed with the kinematically coupled scheme 
with time step At = 1CP 4 (dash-dot line), implicit scheme used by Quaini in [5] with the 
time step A< = 1CP 4 (dashed line) and our scheme with the time step At = 10~ 4 (solid 
line). 
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Figure 9: Example 1: Propagation of the inlet pressure pulse in terms of displacement, 
pressure, and velocity profiles versus time, evaluated at the mid-point of the tube. 
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erence solution and the solutions obtained using At = 5x10 ,10 , 5x10 
and 10~ 4 . Figure 11 shows first-order in time convergence for the velocity. 
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Figure 10: Example 1: Propagation of the pressure wave. 



pressure, and displacement obtained by the kinematically coupled scheme, 
monolithic scheme, and our scheme. The error of our method is notica- 
bly smaller than the error obtained using the classical kinematically coupled 
scheme, and is comparable to the error obtained by the monolithic scheme. 
The values of the convergence rates for pressure, velocity, and displacement, 
calculated using the kinematically coupled schemes, are shown in Table |4j 

5.1.1. Homogeneous Dirichlet vs. absorbing boundary conditions 

We give a short remark related to the impact of the homogeneous Dirichlet 
vs. absorbing boundary conditions. Although absorbing boundary conditions 
for the structure are more realistic in the blood flow application, they will 
only impact the solution near the boundary, except when reflected waves 
form in which case the influence of the boundary conditions is felt every- 
where in the domain. It was rigorously proved in that in the case of 
homogeneous Dirichlet inlet /outlet structure data r\ = 0, a boundary layer 
forms near the inlet or outlet boundaries of the structure to accommodate 
the transition from zero displacement to the displacement dictated by the 
inlet/outlet normal stress flow data. It was proved in [65] that this bound- 
ary layer decays exponentially fast away from the inlet /outlet boundaries. 



Figure 12 depicts a comparison between the displacement obtained using 
absorbing boundary conditions, and the displacement obtained using homo- 
geneous Dirichlet boundary conditions, showing a boundary layer near the 
inlet boundary where the two solutions differ the most. 

It is worth pointing out, however, that absorbing boundary conditions 
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Figure 11: Example 1: Figures show relative errors compared with the kincmatically 
coupled scheme which is first-order accurate in time. Top left: Relative error for fluid 
velocity at t=10 ms. Top right: Relative error for fluid pressure at t=10 ms. Bottom: 
Relative error for displacement at t=10 ms. 



help in reducing reflected waves that will appear when the propagating wave 
reaches the outlet boundary and reflects back. The "optimum" absorbing 
boundary conditions would have to be designed on the basis of Riemann 
Invariants (or characteristic variables) for the hyperbolic problem modeling 
wave propagation in the structure. Those conditions, however, are not always 
easy to calculate, and so approximate Riemann Invariant-based absorbing 
conditions such as (51 )-(52 ) are used. Figure 13 shows the displacement 
in the case of absorbing boundary conditions (51), (52), and homogeneous 
Dirichlet boundary conditions r) = at t = 100 ms and t = 200 ms. Notice 
how the two solutions differ everywhere in the domain, and how the solution 
with absorbing boundary conditions reduces the amplitude and formation of 
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Figure 12: Example 1: Displacement of the structure in the case of absorbing boundary 
conditions (solid line) and homogeneous Dirichlet boundary conditions (dashed line). 
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Figure 13: Example 1: Displacement of the structure in the case of absorbing boundary 
conditions (solid line) and homogeneous Dirichlet boundary conditions (dashed line). 



5.1.2. Example lb. 

This example is an extension of the benchmark problem by Formaggia 
et al. [M] studied in Example 1. The extension concerns inclusion of the 
longitudinal displacement in the model described in Example 1. 

Namely, here we explore how the Koiter shell model ([8])-(|9| looks for the 
coefficients given by those in Example 1. Namely, the radial displacement 
satisfies the same model equation as in Formaggia et al. [M] , while the longi- 
tudinal displacement satisfies equation (fSj) in the Koiter shell model, with the 
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Figure 14: Example lb: Longitudinal displacement (red) and radial displacement (blue) 
for the Koiter shell model in Example 2 calculated with At = 10 -4 . 
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Figure 15: Example lb: Flow rate for the Koiter shell model in Example lb calculated 
with At = 1CT 4 . 
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Figure 16: Example lb: Mean pressure for the Koiter shell model in Example lb calculated 
with At = 1(T 4 . 



corresponding coefficients in accordance with equation (50). More precisely. 



by comparing equations feh and (50) we observe that only the coefficients 



Co, C\ and D\ are different from zero. This implies the following Koiter shell 
model: 
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(54) 
(55) 



Notice that this problem is degenerate in that the operator associated with 
the static equilibrium problem is no longer strictly elliptic. This is due to 
the fact that the coefficient C3 at the second-order derivative with respect to 



z of r] z in equation (54) is equal to zero. Nevertheless, we solve the related 
FSI problem with zero Dirichlet boundary data 77 = 0, and study the flow 



driven by the time-dependent pressure data (53) given in Example 1. The 



values of the coefficients in the Koiter shell model (55)- (54) are equal to 
those in Example 1. Figures [141 [151 and [16] show the displacement, flow 



rate, and pressure, respectively. It is interesting to notice, as is shown in 



Figure 14, that the magnitude of longitudinal displacement is the same as 
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the magnitude of radial displacement. 



5.2. Example 2. 

In this test case we consider the full Koiter shell model ([8])-(|9]). This 
means that all the coefficients Cj and Di are given in terms of the Young's 
modulus E and Poisson ratio <r, and their viscous counterparts E v and a v , 
through the relationships ( Jloj ). Since the terms containing the 4th and 5th 
order derivatives are negligible (it can be shown, using non-dimensional anal- 
ysis, that these terms are much smaller than the remaining terms), we ig- 
nore these terms in the numerical simulation. We present this example as 
a benchmark problem for FSI studies in which the structure is modeled by 
the Koiter shell model (J8|-(10), which captures both radial and longitudinal 
displacement. 

The fluid and structure parameters are given in Tables [I] and |5j The 
model and the coefficients are summarized as follows: 
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where here we take the value of D\ to be equal to 7 from Examples 1 and lb, 
which, from the definition of coefficient D\ above, implies D v = 6R 2, y/h 3 . 
From here we determine C v = D v /a v = 2D V . The corresponding values of 



Parameters 


Values 


Structural viscosity C v (poise cm) 


30 


Structural viscosity D v (poise cm) 


15 



Table 5: Structural viscosity parameters for Example 2. 

the coefficients are given in Table [6] 

This model includes the coupling terms between the longitudinal and 
radial components of the displacement through C 2 7^ and D 2 7^ 0, and 
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Co 


= 4.0133 x 10 5 




= 333.3 


c 2 


= 10 5 


c 3 = io 5 


Do 
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= 10" 2 


D 2 


= 3 


D 3 = 3 



Table 6: Koter shell model coefficients for Example 2. 
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Figure 17: Example 2: Longitudinal displacement rj z , and radial displacement r\ r cal- 
culated on a coarse mesh (solid line) and on a fine mesh (dashed line), obtained with 
At = 10~ 4 . 



the leading-order viscoelastic effects in the radial displacement described by 
Dq 7^ 0. Notice a much smaller value for the coefficient C\ than in Example 1. 
Also notice the large coefficient C 2 that describes the coupling between the 
radial and longitudinal components of the displacement in the Koiter shell 
model. 



Figure [17] shows longitudinal and radial displacement of the structure 
computed on two different meshes, evalued at times t = 2,4,6,8, 10 and 12 
ms. The coarser mesh is twice as fine as the mesh used in Example 1, so 
that the triangularization of the coarser mesh is h°° arse = h p /2, h c ° aTse = 
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Figure 18: Example 2: Flow rate computed on a coarse mesh (solid line), and on a fine 
mesh (dashed line), obtained with At = 10~ 4 . 



h v /2. The fine mesh in this example is twice as fine as the coarse mesh, 
namely h^ ne = h™ arse /2 = h p /A, h{ ine = h c ° arse /2 = h v /A. Note that the 
longitudinal displacement is of the same order of magnitude as the radial 
displacement. Figures 18 and 19 show the corresponding flow rate and mean 
pressure. The appearance of oscillations near the inlet is physical, and is 
associated with a very small value of the coefficient C\ when compared to 
Example 1. An increase in the value of visoelasticity parameters, examined 
in the next example (physiologically reasonable), dampens the oscillations 
observed in this example. 

To study convergence in time we define the reference solution to be the 
one obtained with At = 10~ 6 , and we compute the L 2 -norms of the differ- 
ence in the pressure, velocity and displacement between the reference solution 
and the solutions obtained with At = 10~ 4 ,5 x 10~ 5 , 10~ 5 and 5 x 10~ 6 . A 
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Figure 19: Example 2: Mean pressure computed on a coarse mesh (solid line), and on a 
fine mesh (dashed line), obtained with At — 10~ 4 . 



comparison of the time convergence between our scheme and the kinemati- 
cally coupled scheme is presented in Table [7} Figure 20 shows a comparison 
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Table 7: Example 2: Convergence in time calculated at t = 8 ms. The numbers in the 
parenthesis show the convergence rate for the kinematically coupled scheme presented 
in[S]. 
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between the time convergence of our scheme and the kinematically coupled 
scheme. Similarly as before, our method implemented with (3 = 1 provides 
a gain in accuracy over the classical kinematically coupled scheme, which 
corresponds to (3 = 0. 
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Figure 20: Example 2: Figures show relative errors compared with the kinematically 
coupled scheme which is first-order accurate in time. Top left: Relative error for fluid 
velocity at t=8 ms. Top right: Relative error for fluid pressure at t=8 ms. Bottom: 
Relative error for displacement at t=8 ms. 



5. 3. The Common Carotid Artery ( CCA ) Example 

We conclude this manuscript by showing a simulation of flow and dis- 
placement in the left common carotid artery. A straight segment of the left 
CCA, such as the one shown in Figure 21 , is considered. The geometric pa- 
rameters, such as length, average radius, and wall thickness, are taken from 
the measurements reported in [67J EH EHl El E], while the Youngs modulus 
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Figure 21: A B-mode ultrasound image of CCA [66 representing our computational do- 
main. The black triangle above the CCA is the jugular vein (not a part of the computa- 
tional domain). 



is taken from the measurements reported in [70J. Blood vessels are essen- 
tially incompressible and therefore have the Poisson's ratio of approximately 



0.5 [72]. The Table in Figure 21 shows the values of the corresponding param- 
eters that were used in our simulation. They are within the ranges reported 
in the above-mentioned literature. 



Parameters 


CCA 


Radius H (cm) 


0.3 


Fluid density pf (g/cm 3 ) 


1.055 


Fluid viscosity \x (g/(cm s)) 


0.04 


Wall density p s (g/cm 3 ) 


1.055 


Wall thickness h (cm) 


0.07 


Young's mod. ^(dynes/cm 2 ) 


2 x 10 6 


Poisson's ratio a 


0.5 



Table 8: Geometry, fluid and structure parameters for the common carotid artery example. 



The structural viscosity constants C v and D v are equal to 

C v := 3 x 10 4 dynes/cm 2 • s, D v := C v a. 

This choice of structural viscosity parameters was shown in [45J to resemble 
the viscous moduli of blood vessels. These parameters give rise to the values 
of the Koiter shell coefficients given in Table |9j 
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C = 


1.7022 x 10 6 


Ci = 


846.9 


c 2 


= 3.1 x 10 5 


C 3 = 1.867 x 10 5 


c 4 


= 


D = 


23439.2 


D x = 


9.527 


^2 


= 3500 


D 3 = 2100 


£>4 


= 



Table 9: Koiter shell model coefficients for the common carotid artery example. 



We study blood flow driven by the inlet and outlet pressure data shown 
The morphology of the pressure wave in the left CCA was 



in Figure 22 




0.4 0.6 
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Figure 22: The inlet and outlet pressure data. The average pressure drop is around 0.0673 
mmHg per centimeter. 



obtained from [IS]. The average pressure drop equals around 0.0673 mmHg 
per centimeter, which produces the local Reynolds number of around 1000. 
This is associated with the maximum blood flow velocity of around 100 cm/s, 
which is typical for CCA [731 EH [75], [76] . Indeed, the results of our numerical 
simulation, presented in Figure [23] show the velocity ranging between 22 



cm/s and 97 cm/s, which is within the expected values for the CCA. 

Radial wall displacement has been well examined by many experimental 
studies [7TJ [77J [TSJ, [79] [75]. Maximum radial displacement decreases with 
age, and usually varies between 0.1 mm and 0.38 mm, i.e. between 3%and 



13% of the vessel's radius. Indeed, our simulation, shown in Figure 23 top 
left, indicates maximum radial displacement around 6%, which is well within 
the normal range. 



Figure [23] top right, shows longitudinal displacement computed using our 
thin shell model. We see that the longitudinal displacement in our simu- 
lations lies between min r] z = —0.1 mm and max T) z = 0.05 mm, which 
implies the total longitudinal displacement (defined to be max r] z — min rj z ) 
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Figure 23: From top left, to bottom right: radial displacement, longitudinal displacement, 
pressure, and velocity, calculated at the mid-point of the CCA segment, vs. time. 



of 0.15mm. This is in good agreement with experimental studies obtained us- 
ing B-mode ultrasound speckle tracking method and/or B-mode ultrasound 
velocity vector imaging [801 H El H], which report the total longitudinal dis- 
placement in a healthy CCA between 0.052 mm and 0.302 mm. 

Finally, we observe the captured viscoelastic properties. The viscoelastic 
effect is visible in the stress-strain relationship of the arterial wall, which 
exhibits hysteresis. Figure 24 shows hysteresis between the vessel diameter 
and pressure at the center of the vessel over one cardiac cycle. Hysteresis 
can be quantified by the energy dissipation ratio (EDR), which is a measure 
of the area inside the diameter-pressure loop relative to the measure of areas 
inside and under the loop, i.e., EDR = A 1 /(A 1 + A 2 ) x 100%. Walls with 
higher viscoelasticity have larger area inside the loop, resulting in higher 
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Figure 24: Hysteresis between vessel wall diameter and pressure at the center of the vessel, 
over a cardiac cycle. 

EDR. In our simulations EDR is 8.5%. This is comparable to the results 
in [1H] which show EDR of 7.8% for a young subject. 

6. Conclusions 

In this manuscript we proposed a new thin structure model capturing 
radial and longitudinal displacement of arterial walls, and have designed a 
modification of a loosely coupled partitioned scheme (the kinematically cou- 
pled scheme to numerically simulate the resulting fluid-structure inter- 
action problem between blood flow and arterial walls. The proposed arterial 
wall model is given by the linearly viscoelastic, cylindrical Koiter shell model. 
The fluid and structure are fully coupled using the kinematic and dynamic 
coupling conditions. The new loosely coupled scheme (the kinematically 
coupled /3-scheme), which is lst-order accurate in time, and 2nd-order accu- 
rate in space, is based on a modified Lie splitting. In [6J it is shown that 
this scheme is unconditionally stable for all < /3 < 1. Two test prob- 
lems were presented showing that the kinematically coupled /3-scheme has 
accuracy which is comparable to that of the monolithic scheme by Badia, 
Quaini, and Quarteroni [SI [H] while retaining the main advantages of loosely 
coupled partitioned schemes such as modularity, easy implementation, and 
low computational costs (no sub-iterations between the fluid and structure 
sub-solvers are necessary for convergence). We believe that the main rea- 
son for the increase in accuracy is related to the stronger coupling between 
the leading effect of the fluid load, given by the fluid pressure, and struc- 
ture elastodynamics, which are now strongly coupled in the last step (Step 
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3) of the splitting. It remains to be investigated how does the choice of /3, 
which would provide highest possible accuracy of the kinematically coupled 
(3 scheme, depend on the coefficients in the problem. 

Although the methodology discussed in this manuscript was presented 
for 2D problems, there is nothing is the proposed time-splitting scheme that 
depends on the dimension of the spatial domain. The same methodology 
can be applied to 3D problems, while, of course, the implementation of the 
proposed algorithm in 3D would be much more complex in that case. 

Extensions of the proposed methodology to study FSI between an incom- 
pressible, viscous fluid and a thick structure (arterial walls), modeled using 
equations of 2D/3D elasticity (e.g., St. Venant-Kirchhoff model, used, e.g., 
in (23, |28] to model arterial walls), are under way. 

The research presented in this manuscript provides a first step in our effort 
to capture multi-layered structure of arterial walls and their interaction with 
blood flow. In modeling the intima-media/adventitia complex, the coupling 
between a thin shell (intima) allowing radial and longitudinal displacement, 
and a thick structure (media/adventitia) is important. Development of the 
model presented in this manuscript is a crucial first step. Our preliminary 
results show that the modified kinematically coupled scheme proposed in this 
manuscript is perfect for the numerical solution of such a complex, multi- 
physics FSI problem. Research in this direction in under way. 
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